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We begin by reviewing a technique to approximate the dynamics of stochastic programs — written 
in a stochastic process algebra — by a hybrid system, suitable to capture a mixed discrete/continuous 
evolution. In a nutshell, the discrete dynamics is kept stochastic while the continuous evolution 
is given in terms of ODEs, and the overall technique, therefore, naturally associates a Piecewise 
Deterministic Markov Process with a stochastic program. 

The specific contribution in this work consists in an increase of the flexibility of the translation 
scheme, obtained by allowing a dynamic reconfiguration of the degree of discreteness/continuity of 
the semantics. 

We also discuss the relationships of this approach with other hybrid simulation strategies for 
biochemical systems. 

1 Introduction 

Models in systems biology tend to cluster around two families of mathematical tools: differential equa- 
tions and stochastic processes. Even though, physically speaking, stochastic models have firmer grounds 
[14, 15], their computational analysis is much more costly than that of their differential counterpart. In 
any case, ODE-based descriptions of biological systems are often valuable and provide deep insights. 
Indeed, it is known that, limiting to mass action models, ODE's are an approximation of the average 
of stochastic models, and the differences between the two vanish in the thermodynamic limit [13] (i.e. 
when populations and system's size go to infinity). Recently, there have been many attempts to mix 
these two techniques, at least as far as simulation of biological systems is concerned, resulting in several 
hybrid simulation algorithms [32, 25]. Hybrid dynamical systems have also been a hot topic in the last 
two decades, with much research work spanning across the boundary between computer science and 
engineering control. The best known model among hybrid dynamical systems are hybrid automata [21]. 
Stochastic extensions of such concept are also receiving recently much attention [10], although stochas- 
tic hybrid systems have a somewhat longer tradition [12]. In both cases, most of the interest is in the 
development of automated reasoning tools rather than in simulation. 

It is widely recognized that Computational Systems Biology can highly benefit from modeling 
approaches embodying some stochastic ingredient. A very popular line along which such incorpo- 
ration is realized, is based on the use of stochastic process algebras [27, 11], which are proposed as 
front-end languages to (automatically) generate mathematical models, usually Continuous Time Markov 
Chains(CTMC), see [32]. Recently, such process algebra based languages have also been endowed with 
semantics based on ODE [22], which increase the flexibility of such tools. 

Many proposals of hybrid simulation algorithms for systems of biochemical reactions have been put 
forward [19, 24, 20, 28, 23, 16, 2]. Their salient feature is a description of one part of the system as 
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continuous, while keeping the other discrete and stochastic. The basic idea is to find the best trade off 
between accuracy and computational efficiency (stochastic simulations are much more expensive than 
ODE simulation). 

In this paper we continue a programme which aims to increase even more the flexibility of stochastic 
process algebras by providing them with a very general semantics based on (stochastic) hybrid systems, 
encompassing CTMC and ODE as special cases. Such an approach is motivated not only by the gain 
in flexibility, but also by the goal of exploiting, in a systematic manner, automated reasoning tools to 
provide as much information as possible from a given model. Our stochastic process algebra of choice is 
stochastic Concurrent Constraint Programming (sCCP) [6], an extension of CCP [29] in the stochastic 
setting. In addition to the standard CTMC-based semantics, we have also provided sCCP with an ODE- 
based semantics [4] and with an hybrid automata based semantics. Moreover, hybrid semantics has been 
proposed both with a fixed or user-defined amount of continuously approximated components (see [8, 9]). 

In this paper we extend our work by introducing a semantics based on Stochastic Hybrid Automata, 
thereby guaranteeing the possibility of parameterizing the degree of continuity introduced in the model. 
The approach allows also a dynamic reconfiguration of such degree, in accordance to properties of the 
current state of the system. This allows the description in a formal setting of different hybrid simulation 
strategies, opening the way for their use in the context of process algebra modelling. 

We will start our presentation by introducing, in Section 2, a high level description of the target 
stochastic hybrid systems, suitable to be easily mapped to the well-established formalism of Piecewise 
Deterministic Markov Processes (see supplementary material [1]). The formalism introduced in Sec- 
tion 2, called Transition-Driven Stochastic Hybrid Automata (TDSHA), will act as the intermediate 
layer in the definition of the stochastic hybrid semantics of sCCP. Section 3 briefly introduces the sCCP 
language, while Section 4 presents the mapping from sCCP to TDSHA. Collections of TDSHAs can 
be organized in a lattice, whose definition and basic properties are presented in Section 4.1. Finally, in 
Section 5 we introduce the dynamic reconfiguration mechanism, briefly discussing also how to render, 
in such reconfigurations, partition strategies developed for hybrid simulation algorithms. 

2 Transition-driven Stochastic Hybrid Automata 

We define here a stochastic variant of Transition-Driven Hybrid Automata, introduced in [9] as an inter- 
mediate layer to map sCCP into hybrid automata. The emphasis is on transitions which, as always in 
hybrid automata, can be either discrete (corresponding to jumps) or continuous (representing flows act- 
ing on system's variables). The stochastic variant defined below contains two kind of discrete transitions: 
instantaneous — as in [9] — and stochastic, which happen with an hazard given by a rate function. 

Definition 2.1. A Transition-Driven Stochastic Hybrid Automaton (TDSHA) is a tuple 

ST = (e,X,T£, T£>, 16, init), where: 

• Q is a finite set of control modes. 

• X = {X\,. . . ,X n } is a set of real valued system's variables 1 . 

• %<t is the set of continuous transitions or flows, whose elements x are triples (q, stoich, rate), 
where: q G Q is a mode, stoich is a vector of size |X|, and rate : W 1 — > R is a (sufficiently smooth) 
function. The elements of a triple x are indicated by cmodefr], stoich [t], and ratefr], respectively. 

'Notation: the time derivative of Xj is denoted by Xj, while the value of Xj after a change of mode is indicated by X'- 
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• ID is the set of instantaneous transitions, whose elements 8 are tuples of the form 

(qi , q 2 , priority, guard, reset), where: q\ is the exit-mode, q 2 is the enter-mode, priority : R" — ► R + 
is a weight function used to resolve non-determinism between two or more active transitions. 
Moreover, guard is a quantifier-free first-order formula with free variables in X, representing the 
closed set Gs = {x G R n | guard[x]} in which thew transition is active, and reset is a deterministic 
update of the formX' = /(X). 2 The elements of a tuple 5 are indicated by ei[5], e2[5], priority [S], 
guard[5], and reset[5], respectively. 

• 16 is the set of stochastic transitions, whose elements T7 are tuples of the form 

r\ = (qi,q 2 , guard, reset, rate), where q\, q 2 , guard, and reset are as for transitions in IS), while 
rate : R" — > R + is the rate function giving the hazard of taking transition r\. Such function is 
referred to by rate [17]. 

• init is a point giving the initial state of the system. 

A TDSHA has three types of transitions. Continuous transitions represent flows and, for each T G 
T<£, stoichfr] and rate[r] give the magnitude and the form of the flow of t on each variable XeX, 
respectively (see below). Instantaneous transitions represent actions happening immediately when their 
guard becomes true. Finally, stochastic transitions happen at a specific rate. Both instantaneous and 
stochastic transitions can change system variables according to a specific reset function, depending on 
the variables'value at the point in time at which the jump occurs. 

Remark 2.1. Both priority and rates introduced in Definition 2.1 make TDSHA stochastic. Priorities 
define, at each point, a discrete distribution of a random variable choosing among enabled instantaneous 
transitions. Rates, on the other hand, define a random race in continuous time, giving the delay for the 
next spontaneous jump. 

Product of TDSHA. Given two TDSHA & x = (gi , Xi, Id , TSi , T61 , init x ) and 
^2 = (22, X 2 , T£ 2 , £©2, init 2 ), the product TDSHA 2f = 3\ <g> £? 2 can be defined in a simple way, 
along the path outlined in [9]. Essentially, the discrete states'space of the product automaton is Q\ x Q 2 , 
while transitions from state (qi,qi) are all those issuing from q\ or q 2 . Instantaneous or stochastic 
transitions of 2?\ going from state q\ to state q\ , will go from a state (qi , q 2 ) to (q' l , q 2 ) for each q 2 £ Q 2 . 
Symmetrically for transitions of 

Dynamics of TDSHA. In order to formally define the dynamical evolution of TDSHA, we can map 
them into a well-studied model of Stochastic Hybrid Automata, namely Piecewise Deterministic Markov 
Processes [12]. In this sense, TDSHA are related to communicating PDMP [31], as they can also be 
seen as a compositional formalism to model PDMP. Due to space constraints, we just sketch here an 
informal description of PDMP. The interested reader can find a more formal treatment of PDMP and of 
their relation with TDSHA in the supplementary material [1]. 

Basically, PDMP are stochastic processes whose state space is given by a finite collection of dis- 
crete modes and by a set of real-valued variables. Within each mode, the continuous variables evolve 
following the solution of a set of mode-specific ODE's. While in a mode, variables must stay within the 
allowed region. If they touch the boundary of the allowed region, a forced discrete transition is taken, 
and the system may change mode and/or reset the value of the variables. Moreover, the system is subject 
to the happening discrete stochastic events, governed by an hazard rate that is function of the discrete 



2 Even though there is no real additional difficulty in considering stochastic resets — i.e. in assuming reset to be a transition 
measure — we decided to avoid such move for the sake of simplicity. 



66 



Hybrid semantics of Stochastic Programs with Dynamic Reconfiguration 



mode and of continuous variables. Also stochastic transitions trigger a reset of the state of the system. 
The main points of the mapping from TDSHA to PDMP are the following. 

• Within each discrete mode q G Q, the system follows the solution of a system of ODE, constructed 
combining the effects of the continuous transitions T acting on mode q. Essentially, the ODE for 
variable X, is obtained by adding up the rate of all such t times the i-th component of the vector 
stoichfr]: 

X = ^ stoich[r]rate[T] in mode g G Q. 

t[ |cmode[T]=g 

• Two kinds of discrete jumps are possible: stochastic transitions are fired according to their rate, 
while instantaneous transitions are fired as soon as their guard becomes true. In both cases, the 
state of the system is reset according to the policy specified by reset. Choice among several 
active stochastic or instantaneous transitions is resolved probabilistically according to their rate or 
priority, see Remark 2.1. 

• A trace of the system is therefore a sequence of (random) jumps interleaved by periods of contin- 
uous evolution. 

3 Stochastic Concurrent Constraint Programming 

In this section we briefly present (a simplified version of) stochastic Concurrent Constraint Programming 
(sCCP [3], a stochastic extension of CCP [30]), as it seems to be sufficiently expressive, compact, and 
especially easy to manipulate for our purposes 3 . In the following we just sketch the basic notions and 
the concepts needed in the rest of the paper. More details on the language can be found in [3, 6]. 

Definition 3.1. A sCCP program is a tuple srf = (A, 3>,X,init(X)), where 

1. The initial network of agents A and the set of definitions & are given by the following grammar: 

= | | {C = M} 

% = [ g (x) M (X,X')U( X ) M = n.C\M + M A = M \ A || A 

2. X is the set of variables of the store (with global scope); 

3. init(X) is a predicate on X of the form X = xo, assigning an initial value to store variables. 

In the previous definition, basic actions are guarded updates of (some of the) variables: g(X) is a 
quantifier-free first order formula whose atoms are inequality predicates on variables X and w(X,X') is a 
predicate on X,X' of the form X' = /(X) (X' denotes variables of X after the update), for some function 
f :W —> W. Each such action has a stochastic duration, specified by associating an exponentially 
distributed random variable to actions, whose rate depends on the state of the system through a function 
X :X^M+. 

Example 3.1. We will illustrate the notions introduced in the paper by means of an example coming 
from biological systems. Specifically, we consider a simple model of a (procaryotic) genetic regulatory 
network with a single gene, expressing a protein acting, after dimerization, as a repressor of its own 

3 There are other probabilistic extensions of CCP studied in literature, like [26, 18, 17]. [26] provides CCP with a semantics 
based on discrete time Markov Chains, while in [18, 17] the stochastic ingredient is introduced by extending the store with 
random variables and adding a primitive for sampling. These approaches, however, are not suited for our purposes, as we need 
a model in which events happen probabilistically in continuous-time, as customary in biochemical modeling. 
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production. We assume a cooperative repression: two dimers are required to bind to the promoter region 
of the gene. The sCCP model is given by srf = (A,£t,X,init), where the variables are X = {X p ,X p 2}, 
storing the quantity of the protein p and of its dimer p2 and the components in @ are (* stands for true): 
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[* 




X P + l]t,i • S ene o + [X P 2 > -► *]vx p2 • genei 
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def 


[* 
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X p -2 AX' p2 = X p2 + i]k x x p (x p -l)/2 • dimer + 






[* 




X p + 2AX' p2 = X p2 - l] k _ xXp2 . dimer 



The initial network A is geneo || deg || dimer with initial values of the store variables are given by 

init(Xp,X p2 ) = (X p = 0) A [X P 2 = 0). 

Notice: there is no need to introduce agents for proteins or dimers, as the quantity of these objects 
needs only to be measured by stream variables. The repression mechanism is represented by a gene 
unable of expressing a protein whenever in state gene2- We did not decrement X P2 before entering states 
genei and gene2 as we assume repression mechanism not requiring a binding of the dimer (inhibition 
by bumping). 

Remark 3.1. The pros and cons of using sCCP as a modeling language for biological systems are dis- 
cussed in detail in [6]. Basically, sCCP combines on one side the logical simplicity of process algebras 
and on the other side the computational power of constraints. As a matter of fact, the constraint store 
can be more general than that used in this paper, whereby more complex information (like spatiality) 
can be managed just by a simple programming activity. Further work is needed, however, to export the 
techniques developed here to a more general version of the store. 

All agents definable in sCCP, i.e. all agents C = M£ S>, A are sequential, i.e. they do not contain any 
occurrence of the parallel operator, whose usage is restricted at the upper level of the network. 
sCCP sequential agents can be seen as automata synchronizing on store variables and they can be con- 
veniently represented as labeled graphs, called Reduced Transition Systems (RTS) (see [5]). 

The steps to obtain an object suitable to our subsequent treatment are the following: 

1. Define the collection of all possible states — the derivative set Der(C) — and actions — action(C) — 
of any sequential agent appearing in a sCCP program. 

2. Restrict to sCCP simple programs, i.e. programs without multiple copies of the same agent run- 
ning in parallel at the same time. Formally, it is required that the derivative sets of any two agents 
in parallel in the initial network are disjoint. This is only an apparent restriction, cf. [9] for a more 
detailed discussion. 

3. Introduce the following multi-graph 5 RTS(C) = {S(C),E(C),£): 

• S{C)=Der{C), 

• E{C) = {(exit(7t),enter(7t)) \ % G action(C)}, 

4 In the following, with a slight abuse of notation, we sometimes write C € 2# for C = M e &. 

5 exit(n) ,enter(n) ,guard(n) ,update(n) ,rate(n) give the executing agent, the target agent, the guard, the update and the 
rate of an action n, respectively. 
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• 1(e) = (guard(%) ,update{%) ,rate(7l)), where % is the action defining e G E(C). 
In Figure 1, we show the RTS for the agent geneo, defined in Example 3.1. 
4. Introduce the notion of extended sCCP program 

= (A+,^+,XU{P c \ C£@},init + (XU{P c \ Ce &>})), 

in which a variable Pq for run-time recording the number of parallel copies of each agent C G ^ is 
available, and prove sd + is isomorphic to srf (see [9] for further details). 

Essentially, the last step is a technical trick that simplifies the overall treatment. The variable Pq counts 
the number of copies of C present in parallel within the system at a given point in time. To take into ac- 
count the effects of transitions on agents, we modify updates and rate functions, by increasing/decreasing 
counter Pc relative to actions adding/removing a copy of C. The reason for introducing state variables 
will be apparent in next section. They are required to control a cluster of discrete states (continuously 
approximated) and the real value of a state variable will indicate the "tendency" of the system to be in 
that particular state. 




Figure 1: Reduced Transition Systems for the agent geneo defined in Example 3.1. Each edge is labeled 
by its rate function and by the guard and the update (depicted in the guard — > update notation). 



4 From sCCP to TDSHA 

In this section we define a semantics for sCCP in terms of TDSHA. The basic idea is to partition all pos- 
sible transitions executable by an sCCP agent into two classes: those remaining discrete-stochastic and 
those to be approximated as continuous. Different partitions'schemata correspond to different TDSHA. 
By parameterizing upon such schemata, we will obtain a lattice of different TDSHAs. 
Note that this approach is different from [9], as we do not remove stochasticity . Indeed, Stochastic Hy- 
brid Automata can be seen as an intermediate layer between stochastic programs and (non-stochastic) 
hybrid systems. The reader is referred to [8] for further discussions. 

The mapping proceeds in two steps. First we convert into TDSHAs each sequential component of a 
sCCP program, then all these TDSHAs are combined using the product construction. 

Given a = (A + , @ + , Y, init + (Y)), let C G & + be one of the components of the initial network 
A, and let RTS(C) = (S(C),E(C),£) be its RTS. 

A specific continuous/discrete scheme of approximation is formalized by the choice of a boolean 
vector K G {0, l} m , m = \E(C)\, indexed by edges in£(C): for e G E(C), K[e] = 1 stands for a continuous 
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approximation of the transition, while K[e] = implies that the transition will remain discrete. Let 
E(k,C) = {e£ E(C)\ K[e] = 1} and E(->K,C) = {e £ E(C)\ K[e] = 0}. 

In order to guarantee that the vector field constructed from continuous transitions is sufficiently 
regular, we identify as continuously cipproximable only those actions n such that rate{n) is differentiable 
and rate(n)[X] = whenever guard{n) [X] is false. 6 We call consistent a vector K such that K[e] = 1 
only for edges e that are continuously approximable. In the following, we suppose to work only with 
consistent K. 

At this point we are ready to introduce the basic components of our target TDSHA. 

Discrete Modes. The modes of the TDSHA will be essentially the states S(C) of the RTS(C) . However, 
as continuous transitions cannot change mode, we need to consider as equivalent those states that can be 
reached by a path of continuous edges. Let us denote by ~ K the equivalence relation among states of 
S(C) relating two states if and only if they are connected by a path of continuous edges (i.e. edges in 
E(K,C) of the non-oriented version of RTS(C)). Let S K (C) =S(C)/~ K . For each edge e G E(k,C), we 
define the stoichiometric vector Vy,e as an |Y| -vector, Y = X U {Pc | C G 2$}, such that VY, e [^] = h if and 
only if variable X is updated by transition e according to the formula X' = X + h. 

Example. Consider the gene component of Example 3.1. Its RTS, shown in Figure 1, has three states, 
corresponding to the three components gene , gene 1; and gene 2 , with state variables denoted by Pq, 
P\, and P2, respectively. The RTS has also 6 transitions, indexed by e\,...,e^. Consider the K vector 
equal to (1,0,0,1,1,1): edges 61,64,65, e$ will be approximated as continuous, while the other three 
remain discrete. The relation ~ K has a quotient state space containing two classes: Si = {gene } and 
#2 = {g ene i)g ene 2l- Such a partitioning of the gene's states can be seen as a way to render a slower 
dynamics for the binding/unbinding mechanism of the first repressor, to be compared to a faster one 
relative the second copy of the repressor. 

Continuous flow. The continuous evolution for TDSHA is given by the following set of continuous 
transitions: 

T€ = {([exit(e)],VY, e ,rate(e) \ eeE(K,C)}. 

Stochastic transitions. Stochastic transitions are defined in a very simple way, as guards and rates are 
basically copied from the sCCP edge. The only technicality is the definition of the reset. 

Consider the state counting variables P = {Pc | C 6 They can assume values less than or equal 
to one, as the initial program is simple. Moreover, they range in the whole real-valued interval [0, 1] 
whenever we are in a clustered state [s] collapsing s\ , . . . , Sk of RTS(C). In this case, the state variables 
P Sl , . . . , P Sk must sum exactly to 1 , their value representing the likelihood of state si , . . . , of the cluster 
[s] , respectively. In order to deal with state clusters correctly, we have to ensure that when a state [s] is 
left, all its state variables are set to zero. Moreover, if a discrete transition looping in [s] takes place, then 
the variable of its target state 57 must be set to 1 , while all other variables of [s] are to be reset to 0. To 
enforce this, consider an sCCP edge connecting states s\ and S2, with 

update(e) = X' = /(X) AP^ = P Sl - 1 A/^ = P S2 + 1, 
6 Guards of continuosly approximable n are, in fact, redundant. 
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kuiPi 




P'o = A Pi = 1 

Figure 2: PDMP associated with the gene component of Example 3.1. Transitions approximated contin- 
uously determine a set of ODE, while discrete transitions are stochastic and are depicted here as edges of 
a graph (showing rates, guards and resets). The ODEs are obtained from continuous transitions accord- 
ing to the recipe of Section 2. We chose to display TDSHA in this format as it is similar to the classic 
representation of hybrid automata. 

and define the function fp on P which is 1 on the component corresponding to P S2 and zero elsewhere. 

- /X 

In this way, P' = fp(P) implements the correct updating policy. Let now / combine / and fp: f 

/(X) 
/p(P) 

Putting everything together, we have that the discrete transition associated with e E E(-<K,C) with 
e = (si,s 2 ) is 

([ Sl Us2lguard(e),Y> = f(Y),rate(e)) G 16. 

Instantaneous transitions. At this stage, there is no need to define instantaneous transitions. They 
will be used in Section 5 to deal with dynamic partitioning. 



We can now collect all our considerations into the following definition. 

Definition 4.1. Let srf = (A,@,X,inito) be a simple sCCP program and g/ + = (A + ,@ + ,Y, irat^) be its 
extended version. Let C be a sequential component in parallel in A + , with RTS(C) = (S(C),E(C),£). 
Fix a boolean vector K S {0, l} m , m = \E(C)\. The Transition-Driven Stochastic Hybrid Automaton 
associated with C with respect to K is /7(C, k) = (Q,Y,%<£,'£D,%&,init), where 

• Q = S K (C)=S(C)/~ K ; 

• %€ = {([exit(e)},VY, e ,rate(e) \ e £ E(k,C)}; 

• ID = 0; 

• %e = {([s 1 },[s 2 ],guard(e),Y , = f(Y),rate(e)) \ e = (s h s 2 ) € E(^K,C)}; 

• init = initQ. 



Example. From the previous definition it is easy to generate the TDSHA relative to our running exam- 
ple above, in which K = (1,0,0, 1, 1, 1). Once we have the TDSHA, we can generate the corresponding 
PDMP (see supplementary material [1]), which is shown in Figure 2. 
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P'o = A P'i = 1 



Figure 3: PDMP obtained from the product of TDSHA associated with the three components of the 
sCCP program of Example 3.1. See also the caption of Figure 2. 

Definition 4.1 gives a way to associate a TDSHA with a sequential agent of a sCCP program. In 
order to define the TDSHA for the whole program, we will use the product construction. 

Definition 4.2. Let s& = (A,£F,X,/m'? ) be a simple sCCP program and srf + = (A + , & + ,Y Jnit^) be 
its extended version, with A + = C\ \\ . . . \\ C n . Fix a boolean vector fq for each sequential agent Q. The 
Transition-Driven Hybrid Automaton for the sCCP program with respect to K = (jQ)i=i,-n * s 

3?{A,K) = 5(Ci,/fi)®-'®5(C„,/c n ). 

Example. Consider again the sCCP model of Example 3.1. It has three components: gene, deg and 
dimer, with 6, 1, and 2 edges respectively. We consider three vectors K\ = (1,0,0, 1, 1,2), Kz = (1), and 
K"3 = (1, 1). The product TDSHA of these three components generates the PDMP depicted in Figure 3. 

4.1 Lattice of TDSHA 

Definition 4.2 associates a TDSHA with a sCCP agent for a fixed partition, given by vector K, of the 
transitions into discrete and continuous. Clearly, different choices of K correspond to different TDSHAs, 
with a different degree of approximation of the original sCCP program. The different TDSHAs can be 
arranged into a lattice according to the following pre-order: 

Definition 4.3. Let gfbea sCCP agent, then K\) C ', fc 2 ) if and only if fCi [e] = 1 => K 2 [e] = 1, 

for each transition e G E(A) = E(C\) U . . . UE(C„), with A = C\ \\ . . . \\ C n the initial agent of 

The bottom element of this lattice is obtained for K = 0, while the top element is obtained for K[e] = 1 
if and only if e is continuously approximable. We remind to the reader that transitions not continuously 
approximable must be kept discrete. 

The two "extreme" choices correspond to two particularly important TDSHAs, as shown in the 
following propositions. 

Proposition 4.1. Let stf be a sCCP program. The TDSHA 0) is the CTMC associated with by 

its standard semantics. 

Proposition 4.2. Let g/ be a sCCP program with initial agent A = C\ \\ ... \\ C n . If e is continuously 
approximable for each e £ E(A) = E{C\ ) U . . . U E(C n ), ^=e/, 1) coincides with the system of ODE's 
associated with g/ by its fluid-flow approximation (see [5]). 
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5 Dynamic Partitioning of Transitions 

In the previous sections we have defined a mapping from sCCP to TDSHA fixing the level of discreteness 
and continuity. This choice, however, can be difficult to perform a priori, as one does not know if the 
system will evolve to a state where a different approximation is more accurate. 

This is particularly true when one deals with biological systems. In this case, reactions involving large 
populations of molecular species or having high rates, may be treated as continuous. However, such 
conditions depend on the state of the system and may change during its evolution. Indeed, there has been 
a growing attention on hybrid simulation strategies in systems biology, cf. next Section. 

In order to have dynamic switching, we can extend the discrete modes, introducing states for each 
admissible vector K. New discrete transitions need to be added as well, changing the value of K according 
to some user defined conditions. 

Intuitively, for each transition e the state space must be partitioned in two — possibly empty — regions: 
one where the edge e is treated as a continuous transition and one in which it is discrete. In order to 
define such regions, we consider two predicates: 

1. cont[e] (X), encoding the condition to change edge e from discrete to continuous; 

2. disc[e](X), encoding the condition to change edge e from continuous to discrete. 

An obvious requirement is that the regions identified by cont[e] and disc[e] should be disjoint. In order to 
define such predicates, we will consider a (sufficiently regular, usually continuous) function f e : R" — > R, 
whose sign will discriminate between continuous and discrete regions for edge e. More specifically, we 
define, for a fixed, small e > 0. 

1. cont[e](x) :=/ e (x) > £; 

2. disc[e](x) := f e (x) < -e. 

Using e instead of 0, not only guarantees that the regions in which cont[e] and disc[e\ are true are 
disjoint, but also avoids pathological situations of infinite sequences of instantaneous transitions. See 
supplementary material [1] for further details. 

Now, suppose e is continuous. If the current trajectory enters in a region of the state space in which 
disc[e] becomes true, then we must trigger an instantaneous transition in order to move from K\[e] = 1 
to K"2 [e] = 0. All the variables must remain unchanged. However, in this case it may happen that the new 
relation ~ K2 splits in two the current mode [s] 6 S Kl (C), say [s] Kl = [si] K2 U [S2] k 2 - m this case, we need to 
introduce two instantaneous transitions, one going to [s\] and the other to [S2]. Now, consider the value 
of the state variables of [s], P[ s ] = L.«,g[.s] A,- It can be proved that = 1. Moreover, P^ = +P[ S2 ] 
but, clearly, it is not necessarily the case that the two quantities on the right hand side of the equality 
are equal. This means that the system may "prefer" to move to states in [s\] than to those in [52]. This 
situation is correctly modeled using priorities, i.e. weighting transition to [*,-] by f^.j and re-normalizing 
variables in [s\] and [52] to maintain the property = 1 for each [s] € S K2 (C). 

We now give a formal definition for this construction, following a similar strategy as in Section 4: first 
we construct TDSHA for sequential components, then we apply the product construction to combine the 
local constructions. In order to fix the notation, consider the TDSHA ^{C, K) = (Q, Y, 12) , T6 , init) 
associated with a component C, with respect to a fixed K. With Q K we indicate the set Q K = {([*] K ) \ [s] £ 
Q}. Moreover, T€ K , 13^, and 1& K denote the sets T£, T£>, and T6, respectively, with states in Q 
replaced by the corresponding states (equivalence classes) in Q K . A similar rule applies to init K . 

Definition 5.1. Let si = (A,S>,X,inito) be a simple sCCP program and s/ + = (A + ,3> + ,Y,initQ) be 
its extended version. Let C be a sequential agent in parallel in A + , with RTS(C) = (S(C),E(C),£) and 
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|£"(C)| = m. Moreover, let cont[e],disc[e], e € E{C) be defined as above. The TDSHA with dynamic 
partitioning associated with C is J7(C,cont,disc) = (Q, Y,T£,T2), 16, with: 

1- Q = Uk-g{0,1}'» Qk-, 

2. 1€ = UKe{o,i} m ^K"' 

3. 16 = U Ke { Q1 yn1& K ; 

4. 12) = U Ke {o,i}-^icUl2)o,iUlS)i,o, where 



HD li0 = U[ Sl ] Kl ,[ S2 ] K2 ,P [s2] ,di S c[e]y = g(Y)^ 



e£E(C), 



ffi (e) = 1, K 2 (e) = 0, JCi (V ) = K- 2 (e') /or c/e'^]^ [j 2 ] k 2 / j 

where g assigns value 5^ for s' € [$2] k- 2 , to any other P s , and it is the identity on X. Moreover 

W 

12?0,i = ^[s] Kl ,[s]K 2 ,hcont(e),Y' = Y^j e€E(C), 

K l (e) = 0,K 2 (e) = l,JCi(e') = K 2 {e) fore^e'^j; 

5. init = inti^; 

Definition 5.2. Let srf = (A, @,X,inito) be a simple sCCP program and = (A + ,^ + , Y,/mYj") be its 
extended version, with A + = Ci || . . . || C„. Moreover, fix predicates contj[e],discj[e] for each sequential 
agent Cj of A + , according to Definition 5.1. The Transition-Driven Stochastic Hybrid Automata with 
dynamic partitioning for the sCCP program srf , with respect to (contj,discj)j = i^.. n is 

^A, (contj,discj) y=i,... n ) = ^(Ci,con?i, disci) <8> • ■ ■ <8> 3^C n ,cont n ,disc n ). 

Remark 5.1 (On the fly simulation). In Definition 5.2, the resulting TDHA has a number of modes 
exponential in the number of transitions that sequential agents can perform. This combinatorial explosion 
rules out the possibility of generating all the modes together. However, if we restrict to simulation, this 
is not a real issue, as we need to record only the current mode: the target mode of a transition can be 
generated on the fly as soon as the transition has been taken, given the knowledge of RTS. 



5.1 Hybrid Simulation Strategies 

The hybrid simulation algorithms proposed in literature [19, 24, 20, 28, 23, 16, 2] basically differ in two 
aspects: the kind of continuous dynamics (it can be based on ODE or SDE) and the rules for partitioning 
reactions into continuous and discrete (usually called fast and slow). More specifically, the partitioning 
can be static (done at the beginning of the simulation) or dynamic (i.e. recomputed at run-time). 
Conditions for separating fast and slow reactions are usually twofold: 

1 . the size of species involved in the reaction must all be bigger than a given threshold. Usually, a 
fast reaction j must satisfy a condition like x t > ^T|V;j| for all species i involved in j, where v is 
the stoichiometric matrix. 
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2. the rate function of fast reactions must be reasonably bigger than that of slow reactions. Usually, 
the following constraint is enforced [28]: kj(x)At > A, which ensures that reaction j fires many 
times during the time step At. In [16] a different partition strategy imposes that rates of fast 
reactions are A times faster than the fastest slow reaction, so as to guarantee a separation of time 
scales. 

Dynamical policies sketched above can be easily accounted for in our setting. 
First of all, we need to start from an sCCP model of a biochemical network [6], in which reactions are 
modeled by action capabilities of agents. Then, applying the framework of this paper, we associate a 
TDSHA with such a model, together with a suitable policy for dynamic partitioning of transitions. All 
we have to do is define a function f e for each sCCP transition e, such that f e (x) > when the associated 
reaction can be considered fast and f e (x) < when it is slow. 

As an example, consider a partition strategy based only on the size of populations, like the one 
adopted in [24]. In this case, the function f e for transition e can be the following: 

f e (x) =rmn{xi-K\v[xi,e}\ | v[x h e}^0}, 

where v[-,e] is the stoichiometry of action e, constructed as in Section 4, and A" is a constant (that can be 
tuned for the specific system). Of course, more complex policies can be introduced by suitably modifying 
the functions f e . 

6 Conclusion and Further Directions 

In this paper we provided a specific process algebra, sCCP, with a general semantics based on stochastic 
hybrid systems, parametric with respect to the degree of continuity and discreteness. The different hybrid 
models generated in this way can be arranged in a lattice, and we provided also a way to dynamically 
move within the lattice. This allows to formally describe hybrid simulation algorithms, opening up their 
use as tools to simulate process algebra-based models. Moreover, this approach gives the possibility of 
using other computational analysis methods than simulation, like reachability computations or model 
checking. 

An interesting problem is how to extend such machinery to other process algebras. First steps have been 
done to deal with stochastic 7T-calculus [7] , however the peculiarities of each language present specific 
difficulties to be solved. 

The formal treatment developed in the paper, in particular the lattice of TDSHA defined in Sec- 
tion 4.1, can also provide an interesting theoretical framework to study the quality of the approximation 
and the error introduced. In particular, the mature theory of PDMP [12] can provide interesting tools in 
this direction. 

Another issue we are investigating regards the relationships between discreteness and stochasticity. 
In particular, we are interested in understanding whether the stochastic ingredient of the dynamics can 
be dropped in favor of a pure discrete evolution, and at what price [8]. Motivations for this reside in the 
fact that non-stochastic hybrid systems have a much wider and more efficient set of automated reasoning 
tools available. 

We conclude with a more basic (perhaps philosophical) question: given that a mix of continuous 
and discrete simulation strategy is the choice, is there a way — other than minimization of computational 
complexity — to determine which parts of the systems can/may be simulated discretely/continuously? 
We feel that physical consideration must be taken into account for addressing this issue and that these 
are probably outside our reach. However we wish to contribute the "computer scientist point of view": 
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the level of discreteness/continuity can be established on the ground of a formal specification of the 
properties to verify/simulate and should guarantee the minimum of computational resources necessary 
to this task. 
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